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. In the present paper, we shall study the 4-dimcnsional Z-2 lattice gauge model with a random gauge 

o 

phase transition. We numerically obtain a phase boundary curve in the (p — T)-plane where T is 
the "temperature" measured in unit of J/ks- This model plays an important role in estimating the 



coupling; the random-plaquette gauge model(RPGM). The random gauge coupling at each plaquette 
takes the value J with the probability 1 — p and — J with p. This model exhibits a confinement- Higgs 



accuracy threshold of a quantum memory of a toric code. In this paper, we are mainly interested 
in its "self-duality" aspect, and the relationship with the random-bond Ising model(RBIM) in 2- 
dimensions. The "self-duality" argument can be applied both for RPGM and RBIM, giving the 
same duality equations, hence predicting the same phase boundary. The phase boundary curve 
obtained by our numerical simulation almost coincides with this predicted phase boundary at the 

5_i ■ 

high-temperature region. The phase transition is of first order for relatively small values of p < 0.08, 
but becomes of second order for larger p. The value of p at the intersection of the phase boundary 
curve and the Nishimori line is regarded as the accuracy threshold of errors in a toric quantum 
memory. It is estimated as p = 0.110 ±0.002, which is very close to the value conjectured by Takeda 
and Nishimori through the "self-duality" argument. 
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1 Introduction 



The duality transformations provide us with important informations for certain class of statistical- 
and field-theoretical models. The most famous one is the Kramers- Wannier self-duality relation for 
the 2-dimensional(2D) Ising spin modelpp, which predicts the exact value of the critical temperature. 

Its gauge- model counterpart was studied by Wegner[3Ej, and he showed that the 4-dimensional 
Ising(Z2) lattice gauge theory is also self-dual. Its self-duality condition for the critical coupling con- 
stant is equivalent to that of the 2D Ising spin model, so the value of the critical coupling(temperature) 
coincides with that of the Ising model, though the orders of the phase transitions may be (and ac- 
tually are) different in these two systems @|. 

Recently, a duality transformation is applied for a random spin system of spin glass, the random- 
bond Ising model(RBIM) in two dimensions This model contains a new parameter p, which 
controls the random quenched variables, i.e., the rate of nearest-neighbor spin-coupling with "wrong 
sign" . Thus the critical temperature becomes a function T(p) of p. The "self-duality condition" for 
this model is proposed in order to locate the multicritical point 6 . The multicritical point is the 
intersection of the critical temperature T(p) and the Nishimori line of spin glass. (More detailed 
discussion on this point will be given in Sec. 2.) The duality transformation for the 2D RBIM is exact, 
but contrary to the nonrandom case (j> = 0), the proposed "self-duality condition" does not assure 
us that the singular point of the free energy is located at the "self-dual point" . However, once this 
"self-duality condition" is accepted, it allows us to make a conjecture on the location of T(p). After 
that, this conjecture for RBIM has been verified by numerical simulation^ for high-temperature 
region. 

The 4-dimensional (4D) Z 2 random-plaquette gauge model(RPGM) plays an important role in 
the theory of quantum memory. Here the randomness p is the probability that the gauge coupling for 
each plaquette takes the "wrong-sign". (See Sect. 2 for more details.) This model is used to predict 
the accuracy threshold of a 3-dimensional(3D) toric quantum code[Sl El- Actually, the accuracy 
threshold p c of the errors of a 3D toric quantum memory is determined by the multicritical point of 
the 4D RPGM^J. The accuracy threshold of a 3D toric code p c is expected to be higher than that 
of a 2D toric code, p c = 0.033^21 • Thus it is quite important to obtain the phase boundary of the 
4D RPGM. 

Takeda and Nishimori ^2 EI] applied the duality transformation for the 4D Z 2 RPGM, and 
assumed the "self-duality condition" in order to locate the multicritical point as p c = 0.110028.... 
The "self-duality condition" in the 4D Z 2 RPGM is the same with that of the 2D RBIM, and so 
the phase boundaries T(p) of these two models are the same if the "self-duality condition" correctly 
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predicts the phase boundary. 

In this paper, we shall study the phase structure of the 4D RPGM by numerical simulation, 
in particular its confinement-Higgs phase transition. The result of the phase boundary curve shall 
be used to determine the accuracy threshold of a 3D toric code as well as to judge whether the 
"self-duality condition" proposed to determine T{p) in Ref.^JE] is valid. We calculate the phase 
boundary curve, i.e., the critical temperature T{jp), in the p — T plane. The curve T(p) starts from 
the critical point of the nonrandom gauge model with the uniform coupling constant at p — and 
decreases as the randomness p increases. Our result T(p) is plotted in Fig.l together with the result 
of "self-duality condition" of Ref. jlll IT3"] [See Eq.lJUJ]. The order of the phase transition is of first 
order for p < 0.08, but it becomes of second order at larger values of p. We recall that the similar 
behavior in the order of phase transition is observed in the 3D RPGM, in which the phase transition 
changes from the second order transitions to the higher order ones^2j- Our curve T(j>) of the 4D 
RPGM coincides with the predicted phase boundary by the "self-duality conditions" of the both 
models mentioned above. This result verifies the conjecture by Takeda and Nishimori.il, 13; for the 
multicritical point of the 4D RPGM. 

This paper is organized as follows. In Sec. 2, we briefly review the duality transformation of the 
random models and show how the "self-duality condition" determines the critical curve in the p — T 
plane. In Sec. 3, we report our numerical calculations presenting the specific heat, the expectation 
values of the Wilson loops, and the fluctuations of the internal energy of the 4D RPGM, to determine 
the critical curve. Section 4 is devoted for conclusion. 

2 Duality in the 4D random-plaquette gauge model 

In this section we shall briefly review the duality transformation of the 4D RPGM following Ref. 11 . 
We consider the 4-dimensional hypercubic lattice, and put a Zi gauge variable U Xil — ±1 on each 
link (xfi), where x denotes the site and fi — 1,2,3,4 is the index for four positive directions. The 
energy H of the model is given by 



where J(> 0) is the strength of the (inverse) gauge coupling, J\ p U is the product of the four gauge 
variables U x ^ on the four links surrounding the plaquette P. Tp is the random variable for each 
plaquette P taking rp = 1 with the probability 1 — p and the "wrong-sign" Tp = — 1 with the 
probability p (G [0, 1]). 

We employ the replica technique for taking the ensemble averages over Tp. For the n-replica 
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system, the averaged partition function Z n is given as 

n n ^ 

Z"(r)=(lJ J2 )exp(-/3^i?(C/ Q ;r)), ^e- , (2) 

a=l£/«,=±l a=l 

where a = 1, • • • , n is the replica index, and (0(r)) en s denotes an ensemble average over {rp}. The 
free energy _F is obtained by taking the limit n — > as usual, 

F = -IlinxffclV (3) 



(3 n->o y n 

Let us introduce the Boltzmann weight for a single plaquctte P, k+ = cxp(K) for the fluxless 
configurations Yip U — 1 and the weight k_ = cxp(— K) for the fluxfull ones, where K = (3 J = 1/T. 
As we are considering n replicas, it is useful to introduce the generalized Boltzmann weight for a 
single plaquette, i.e., for the configurations with fluxless n — k replicas and fluxfull k replicas, we 
assign the weight x k as|?H[TT] 

X k = {l- P )Kl- k K k _+pK%K n _- k . (4) 

Then the averaged partition function is expressed in terms of these weights, x%, ■ ■ ■ , x n , 

Z n = Z n (x ,Xi,- ■ ■ ,x n ). (5) 

To make a duality transformation, we introduce the dual Boltzmann weights x* k {k = 0, • • • , n) 
by the following discrete Fourier transformations; 

1 



(K+ + K-) n ~ 2m ( K+ - K_) 2m 



2 n/2 
f 



*2m+l ^(l-2p)( K+ + K -T- 2m - 1 (K + -K-f m+1 . (6) 

Then the following duality relation can be derived [TT]. 

Zn(%0: -El: ' ' ' 7 %n) — Z n {x^ , X-y , ■ ■ ■ , X n ) , (7) 

up to an irrelevant overall constant. 

In the standard (nonrandom) Zi gauge model with p = 0, it is known that a confinement-Higgs 
phase transition takes place at K = Kc where e~ 2Kc = \/2 — 1. This critical value is obtained by 
imposing the self-duality condition xo = Xq, with which the other n conditions x k = x^. (k = 1, • • • , n) 
also hold (automatically) at p = 0. Then it is expected that the phase boundary of the confinement- 
Higgs phase transition evolves starting at (p = 0, T = 1/Kc) into the region < p. It is very 
interesting to see if the "self-duality condition" , 

x = Xq (8) 
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predicts not only the location of the multicritical point but also that of the whole phase boundary. 
Here we note that one cannot impose n self-duality conditions Xk = x\ (k = 1, • • • , n) simultaneously 
for < p; these equations are overcomplete and have no solutions in contrast with the case of p = 0. 
The arguments of the original and transformed partition functions are not equal even at the "self- 
dual point" . Therefore it is not necessarily assured that the "self-duality condition" JH|l determines 
the singular point of Z n . 

Explicitly in the limit n — > 0, Eq.JSJ) reduces to the following equation, 

which determines the phase boundary curve T(p) in the p — T plane. The multicritical point is 
defined as the intersection point of the two curves, the phase boundary T(p) and the Nishimori line 
defined by 

exp(-21T) = (10) 
1-p 

The value p — pc at the multicritical point is regarded as the accuracy threshold for the error rate of a 
quantum memory of the 3D toric code^JEl- By using Eq.Q, Takeda and Nishimori 11 determined 
Pc as pc = 0.110028..., which is considerably larger than the accuracy threshold p c = 0.033 for a 
2D toric code^2> as ^ is naturally expected. 

Here we comment on the 2D RBIM. In Ref.jHl E], the duality transformation has been ap- 
plied for the 2D RBIM. By imposing the same "self-duality condition" as Eq.(JH|, one obtains just 
Eq.(|5J), so these two models are predicted to have the same phase boundaries T(p). The numerical 
simulation|141 17| of the RBIM gives the phase boundary that almost coincides with Eq.© in the 
region of the p — T plane above the Nishimori line (the high-T region). 1 As explained above, the 
"self-duality condition" for the random systems is just a conjecture, so it is interesting to see if it 
is satisfied also in the random gauge systems. The results of the 4D RPGM will be reported in the 
following section. 



3 Numerical study of the 4D Z 2 RPGM 

In this section, we shall show our results of numerical simulation for the phase structure of the 
4D Z2 RPGM, particularly whether the phase boundary of the confincment-Higgs phase transition 
coincides with the predicted curve . As explained in the introduction, there exists a confinement - 
Higgs phase transition in the nonrandom gauge theory with p = 0, and the critical coupling is given 

1 To obtain definite results by the numerical studies for the low—T region is rather difficult and it requires a very 
large number of random samples. 
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as Kc = —\ ln(\/2 — 1) by the self-duality nature of the model. We expect that the phase boundary 
curve evolves from the point (p = 0,T = 1/Kc) downward as p increases in the p — T plane since 
the inclusion of random couplings tends to put the system in a disordered phase. 

In our simulation, we first generate Tp over the lattice randomly to prepare a sample. Then 
we perform Monte Carlo simulation of this sample with Metropolis algorithm. After repeating this 
procedure, we get the results of a set of samples. Finally we average these results over the samples. 
We calculated the following quantities; 

1. Internal energy and specific heat 

2. Expectation values of the Wilson loop and their deviations from the perimeter/area law 

3. Fluctuations of internal energy and specific heat among samples 

The internal energy and the specific heat are useful to determine the phase boundary and the order 
of the phase transition in the high-T region. As we explained before, the high-T region is the region 
above the Nishimori line in the p — T plane. The second quantity, the Wilson loop, is an order 
parameter of the gauge theorv[TS|: it obeys the area law in the confinement phase, while it obeys 
the perimeter law in the deconfinement phase. In the present case, we use it to locate the phase 
boundary close to the Nishimori line in the high-T region. The third quantity is used to identify 
the multicritical point, the intersection of the phase boundary curve and the Nishimori line. As it 
can be proved exactly, the specific heat (and the internal energy) shows no singular behavior on 
the Nishimori line |16l IT2*] . Then it is rather difficult to identify the phase transition point near the 
Nishimori line. We use all the above three quantities to identify the transition points. 

In Fig. 2, we plot the internal energy (E) and the specific heat C per site at p = 0, 0.02, 0.04, 0.06, 0.08 
and 0.10. It seems obvious that the phase transition is of first order for smaller value of p, p < 0.08, 
whereas it becomes of second order for larger p. The critical points T(p) in Fig.l in the high-T re- 
gion are determined from these peaks of C. They are fairly in good agreements with the self-duality 
prediction ©, but there still exist small differences. We think that they are due to the finite-size 
effect. To check this point, we study the size-dependence of the specific heat. Fig. 3 shows C at 
p = 0.08 for the lattice sizes A 4 with N = 6, 12, 16. We observe that the peak of C becomes sharper 
and higher for larger lattices, which verifies the second-order phase transition at p = 0.08. The value 
of T at the peak decreases gradually as N increases. In Fig. 4 the location of T at the peak of C is 
plotted versus the inverse of the linear size N of the lattice. As N increases, T(p) gets closer to the 
predicted value T(p = 0.08) = 1.552... by the self-duality condition (J5J). 

Next let us see the expectation value of the Wilson loop (for the Z-i gauge theory it was first 
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introduced by Wegner 0), W(C), for a close loop C on the lattice, which is defined as follows; 
W{C) = (W(C, T )) cns , 

W(C,t) = m2(U u ^) e M-PH(u^))/z(T), 

x ^ U ml _ t c 

z m = nE°p(-Wi T ))' ( n ) 

In the confinement phase of the gauge theory, W(C) obeys what is called the area law, 

W{C) oc e- aA{c \ (12) 

where a is a constant named string tension and A(C) is the minimum value of the area of a membrane 
that covers the closed loop C. On the other hand, in the deconfincmcnt phase, W(C) obeys the 
perimeter law, 

W{C) oc e- 7P(c) , (13) 

where 7 is another constant and P(C) is the perimeter of C. 

We calculate W(C) for various p's with a fixed K = 1/T and fit the data by the area and 
perimeter laws. In Fig. 5, the typical results near the multicritical point are shown. These results 
clearly show that W(C) changes its behavior from the area law to the perimeter law as p increases. 
In Fig. 6, we present the results of \ 2 fittings for K = 0.9 and 1.0. We can (roughly) estimate the 
phase transition point by using these results. 

Finally, we study the fluctuations of (E) and C over the samples. In the previous studies on 
the 2D RBIM and the 3D RPGM, it was observed that the fluctuations of these quantities indicate 
signals of the phase transitional El- ln Fig. 7, we present these fluctuations for K = 0.9 and 
1.0. The fluctuation of C stays almost constant at large p (i.e., in the confinement phase), and it 
increases considerably as p decreases and passes a certain critical value. The critical point of p is 
very close to the phase transition point estimated by the other numerical calculations given above. 
More precisely, we estimate the values of p at the criticality for K — 1.0 as p — 0.110 ± 0.002. 
Since this point is very close to the expected multicritical point as seen in Fig.l, we think that it is 
reasonable to use this value as the estimation of p c at the multicritical point. 

4 Conclusion 

In this paper we studied the 4D RPGM numerically. All the calculations consistently indicate the 
existence of the confinement-Higgs phase transition continuing from the nonrandom case p = to 
the random case < p. The critical curve T(p) is determined synthetically from (E), C, W(C), 
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and the fluctuations of (E) and C over the samples. We summarize the results in Fig.l, which 
are in good agreement with the predicted value by the "self-duality condition" 0. In particular, 
our estimation of the multicritical point is p c = 0.110 ± 0.002, while the "self-duality condition" 
predicts p — 0.110228... We regard the small but non-negligible differences between the numerical 
result of T(p) and © as the finite-size effect. Our studies strongly indicate the correctness of the 
conjecture by the "self-duality" © not only for the spin glass model but also for the random gauge 
model in the high-T region. Numerical simulations of the random models at low-T region require 
considerably many samples to obtain reliable ensemble averages so we did not present T(p) at 
the low-T region below the Nishimori line in Fig.l. It is reserved for the future problem. 
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Phase Diagram 




Figure 1: Phase diagram of the 4D RPGM and the 2D RBIM in the p(randomness)-T(temperature) 
plane. The diamond symbols show the transition points of the 4D RPGM obtained by our numerical 
study. The cross symbols show the transition points of the 2D RBIM calculated in Ref. ,14 . The 
triangle symbol shows the transition point of the 2D RBIM on the Nishimori line calculated in 
Ref. [7] . The dashed curve is the phase boundary obtained by the "self-duality condition" @ , which 
predicts the same curve for the 4D RPGM and the 2D RBIM. This curve separate the confinement 
phase on the right-hand side (larger p) and the Higgs(deconfinement) phase on the left-hand side 
(smaller p). The solid curve is the Nishimori line (|10fl . which is also identical for the two models. 
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Figure 2: Specific heat C and internal energy (E) per site vs 1/T in the high-T region for typical 
samples (from p = up to p = 0.1). For p < 0.08, C exhibits a double-peak structure, and (E) 
exhibits hysteresis, both indicating that the phase transition is of first order. The lattice size is 
16 4 for p = 0.0 and 12 4 for the other p's, and the typical sweep number is 3 x 10 4 for p — and 
5 x 10 4 (p = 2,4,6%),10 5 (p = 8%), 2 x 10 6 (p = 10%). 
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Specific Heat (error rate = 8% ) 
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Figure 3: Lattice size dependence of the specific heat for p — 0.08. The number of used samples is 
~ 10 2 . 
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Figure 4: Lattice size dependence of the transition point for p = 0.08. The bold line indicates the 
predicted value T = 1.552... by the "self-duality condition" 
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p = 0.105 p = 0.13 

Figure 5: Wilson loops W(C) for typical samples in the region near the Nishimori line. The 
lower(upper) horizontal axis P denotes the perimeter of the loop C for the perimeter-law fitting, and 
the upper(lower) horizontal axis A denotes the area of C for the area-law fitting. The data show the 
perimeter(area)-law behavior for p = 0.105(p = 0.13) in the left(right) column at each temperature. 
The lattice size is 10 4 , and the typical number of sweeps is 2 x 10 4 . 




0.105 0.11 0.115 0.12 0.125 0.13 °- 105 °- 115 ■ 12 °- 125 ■ 1 3 



(a) T- 1 = 0.9 (b) T- 1 = 1.0 

Figure 6: Deviation x 2 over 500 samples from the perimeter or area law vs p. W(C) is normalized 
as — ln[W(C5)] = 1, where C5 is the contour of a square whose area is A(C$) = 5 2 . The arrows 
show the transition points estimated from the specific heat and internal energy (see Fig. 7). 
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Figure 7: Specific heat and internal energy in the region near the Nishimori line. The histograms 
show their fluctuations among 200 samples. The arrows exhibit the phase transition points at which 
the magnitudes of fluctuations change their rate of variation with respect to p. The lattice size is 
10 4 , and the typical number of sweeps is 2 x 10 5 . 
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